Computation of dendritic microstructures using a level set method 
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We compute time-dependent solutions of the sharp-interface model of dendritic solidification in two 
dimensions by using a level set method. The steady-state results are in agreement with solvability 
theory. Solutions obtained from the level set algorithm are compared with dendritic growth simula- 
tions performed using a phase-field model and the two methods are found to give equivalent results. 
Furthermore, we perform simulations with unequal diffusivities in the solid and liquid phases and 
find reasonable agreement with the available theory. 
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Various numerical approaches |l]-[§] have been devel- 
oped to solve the difficult moving boundary problem that 
governs the growth of dendrites |i|-f7| . Unfortunately, the 
direct solution of the time-dependent Stefan problem is 
troublesome and usually requires front tracking and lat- 
tice deformation in order to contain the moving solid- 
liquid interface, which is often very complicated topolog- 
ically. In general, the methods developed to tackle the 
free-boundary problem have difficulty in handling topol- 
ogy changes, such as the merging and breaking of sur- 
faces, and are usually not easily extendible to higher di- 
mensions. 

In order to avoid the difficulties associated with track- 
ing a sharp interface, the phase-field model of solidifica- 
tion has been developed and is currently the most popu- 
lar technique for simulating dendritic growth. The phase- 
field model avoids the computational difficulties associ- 
ated with front tracking by introducing an auxiliary or- 
der parameter, or phase-held, ip(r,t) that couples to the 
evolution of the thermal field. The dynamics of ip(r,t) 
are designed to follow the evolving solidification front 
|p|-|rT|, which is defined by the zero level set ip(r,t) = 0. 
Because the interface is never explicitly tracked, compli- 
cated topology changes are handled easily. Furthermore, 
the extension of the phase-field model to higher dimen- 
sions is straightforward. 

Although phase-held models have been very useful in 
studying solidification patterns, there are still some lim- 
itations in this approach. The proper use of these mod- 
els requires that an asymptotic analysis be performed in 
order to obtain a mapping between the parameters of 
the phase-field equations and the sharp-interface equa- 
tions [p^|-|l4|. The asymptotics involve expanding the 
phase-field equations in some small parameter propor- 
tional to the interface width, W, and as a result, the 
phase-field model only reproduces the dynamics of the 
sharp-interface equations in the limit where the expan- 
sion parameter is sufficiently small. Computationally, the 
grid spacing must be small enough to resolve the interfa- 



cial region, which is on the order of W. This restriction 
is generally not a problem for the symmetric model of 
solidification (where the diffusivities in the solid and liq- 
uid phases are assumed to be the same) because it is 
possible to have W on the order of the capillary length 
flI4| ]. However, phase-field asymptotics for unequal diffu- 
sivities can be problematic |l5[ ; correction terms that are 
inconsistent with the sharp-interface equations are gen- 
erated and non-monotonic behavior is required in the in- 
terfacial region, which requires extra grid resolution and 
hence slower computational performance. The general- 
ization of the phase-field approach to handle discontinu- 
ous material properties requires a better understanding 
of the mapping between the phase-field model and the 
sharp-interface formulation in order to avoid problems 
with properly resolving the interface. 

The level set method is a computational approach that 
has the capability of avoiding the above mentioned limi- 
tations of front tracking methods and phase-field models. 
This method, first introduced by Osher and Sethian jl6) , 
is conceptually similar to a phase-field model in that the 
solid-liquid interface is represented as the zero contour of 
a level set function, </>(r, i), which has its own equation 
of motion. The movement of the interface is taken care 
of implicitly through an advection equation for <fi(r,t). 
Thus, topology changes and the extension of the method 
to higher dimensions can be handled in a straightforward 
manner. Unlike the phase-field model, there is no arbi- 
trary interface width introduced in the level set method; 
the sharp-interface equations can be solved directly and, 
as a result, no asymptotics are required. Discontinuous 
material properties can also be dealt with in a simple 
manner. 

The level set method has been applied to several prob- 
lems involving moving boundaries |17|-[l9||, including so- 
lidification. Prior work on dendritic growth includes an 
application of the method to a boundary integral formu- 
lation 1 20 1 as well as the direct solution of the sharp- 
interface equations [|T). While these simulations have 
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reproduced the qualitative features of dendrites, as well 
as some quantitatively accurate solutions to exactly solu- 
ble problems, some of the simulations of anisotropic den- 
dritic growth were not necessarily converged |2lf| . Fur- 
thermore, the results were not compared with theoretical 
predictions of dendritic growth. 

In this paper, we demonstrate that the level set method 
can be used to solve the free-boundary problem for solidi- 
fication to calculate quantitatively accurate solutions for 
dendritic growth. We present results from simulations 
in two dimensions and show that the solutions converge 
to the steady-state predicted by microscopic solvability 
theory. Time-dependent results are also compared with 
calculations using a phase-field model and good agree- 
ment is found for all times. Furthermore, wc perform 
simulations with unequal diffusivities (a case which is 
not yet possible with phase-field models) and find that 
the prediction of Barbieri and Langer p2| provides a fair 
quantitative fit to our results. 

The solidification of a pure substance is described by 
a free-boundary problem for the temperature in the solid 
and liquid phases, and the position of the interface be- 
tween them: 

d t u = DV 2 u (1) 

V n = (Dd n u) SoHd - (Dd n u) Liquid (2) 
Ui = -d{6)K - P{6)V n (3) 

The temperature T has been rescaled as a dimensionless 
thermal field u = (T — T m )/(L/C p ), where T m , L, and 
C p represent the melting temperature, the latent heat 
of fusion, and the specific heat at constant pressure, re- 
spectively. The thermal diffusivity, D, can be different 
in the solid and liquid phases. Eq. 2 describes energy 
conservation at the solid-liquid interface, where V n is the 
local outward normal interface velocity and d n refers to 
the outward normal derivative at the interface. Finally, 
Eq. 3 is known as the Gibbs-Thomson condition and de- 
scribes the deviation of the interface temperature, Ui, 
from equilibrium due to the local curvature, k, and in- 
terface kinetics. d(9) = -f(6)T m C p / L 2 is the anisotropic 
capillary length, proportional to the surface tension j(9), 
and (3(9) is the anisotropic kinetic coefficient. Here we as- 
sume that (3(9) — and that the capillary length has the 
form d(9) = d (l — 15ecos46'), where e is the anisotropy 
strength and 9 is the angle between the local normal vec- 
tor at the interface, n, and the x-axis. 

We solve the above free-boundary problem by using a 
level set algorithm, which involves the following steps: 
i) advancing the interface, ii) reinitializing the level set 
function to be a signed distance function, and iii) solving 
for the new thermal field. The general level set method is 
described below. We wish to note that in our simulations 
we implement a localized level set method, described in 
detail in Ref. p3], in which calculations of (j> are per- 
formed only in a narrow region around the interface. We 
have not yet made an attempt to make our algorithm 



more computationally efficient by using adaptive mesh 
refinement. 

i) Advancing the interface. The level set function is de- 
fined as the signed normal distance from the solid-liquid 
interface such that <f> is positive in the liquid phase, neg- 
ative in the solid phase, and zero at the interface. (f> 
satisfies the pure advection equation 

^+F\Vcj ) \=Q (4) 

Integrating Eq. 4 for one timestep results in moving the 
contours of <p along the directions normal to the interface 
according to the velocity field F, which varies in space. F 
is constructed to be an extension of the interface velocity, 
V n , such that F — V n for points on the interface and the 
lines of constant F are normal to the interface. Thus, 
advecting </> according to Eq. ^ moves the front with the 
correct velocity. 

Rather than using a partial differential equation to 
generate F (as in Refs. |2^j2^]), we construct F in the 
following manner: cf) represents the normal distance from 
the solidification front, so the value of cf> at each grid- 
point on the computational lattice can be used to locate 
a particular point on the interface. If x g is the location 
of the gridpoint, the associated point on the interface is 
at Xi — Xg — (j)ri, where the normal vector n = V</>/|V0|. 
The temperature at xi is then calculated by using Eq. 3; 
9 is easily found from ri 7 and the curvature, k — V • n, is 
interpolated at Xj from values of k at neighboring grid- 
points, n and k are calculated using standard, centered 
finite difference approximations to the partial derivatives 
of <f). Next, values of u are interpolated in both the liq- 
uid and solid phases, a distance Ax (the size of the grid 
spacing) away from Xi along the normal direction. These 
two interpolated temperatures are used along with Ui to 
approximate the difference in the normal derivative of u 
at Xi and thus find V n (Eq. 3). Because Xi and x g lie in 
the same line normal to the interface, the value of F at 
x g is simply V n . The field F can be determined at all 
gridpoints in this way. 

After F is known, the interface can be advanced one 
timestep. For stability, we discretize Eq. fusing a 5th or- 
der WENO (weighted essentially non-oscillatory) scheme 
in space and a 3rd order Runge-Kutta scheme in time 
p4j . However, the overall accuracy of our algorithm is 
second order in space and first order in time. 

ii) Reinitialization. After solving Eq. ^ for one 
timestep, the level set function will no longer be equal 
to the distance away from the interface. It is necessary 
to reinitialize <j) to be a signed distance function. This 
step is accomplished by solving 

^+<?(0)[|V0|-1] = O (5) 

to steady state. S((j>) takes on the value +1 in the liquid 
phase, —1 in the solid phase, and is zero at the interface. 
We typically iterate Eq. M three times in order to obtain 
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an accurate distance function. Like Eq. |4|, this equation 
is discretized using a 5th order WENO scheme in space 
and a 3rd order Runge-Kutta scheme in time Q . 

hi) Solving for the new thermal field. The thermal 
held is updated by solving Eq. 1 using a modified Crank- 
Nicolson scheme. Different diffusivities in the two phases 
can be taken into account by simply noting the sign of 
the level set function and using the appropriate diffusion 
coefficient in the finite difference stencil. Special care has 
to be taken for gridpoints near the interface. If \<fi\ < Ax, 
the level set function is used to determine whether the 
front intersects the stencil and, if so, interpolate where 
the interface crosses the stencil. The stencil is then mod- 
ified to take into account the location of the interface and 
the Gibbs-Thomson condition. 

We compute four-fold symmetric dendrites iaaLxL 
square box using the procedure described above. Solidifi- 
cation is initiated by a small quarter disk of radius R a in 
the lower left-hand cor ner of the box. The initial level set 
function is (j>(x, y) = y^x 2 + y 2 — R Q , where x and y are 
the usual Cartesian coordinates. The initial temperature 
is u = in the solid and decays exponentially away from 
the interface to u = — A as x — > oo, where the far-held 
undercooling is A = (T m — T 00 )/(I//C P ) and is the 
temperature far ahead of the solidihcation front in the 
liquid. 

Eqs. 1-3 have been studied extensively to determine the 
steady state features of dendritic growth . According 
to microscopic solvability theory, these equations admit 
a family of discrete solutions. Only the fastest growing 
of this set of solutions is stable. This solution is the dy- 
namically selected "operating state" for the dendrite and 
corresponds to a unique tip shape and velocity. Recent 
calculations of dendritic growth using phase-held mod- 
els have been found to be in good agreement with the 
predictions of microscopic solvability theory Jl4|,|25| ■ We 
observe similar agreement with the use of the level set 
algorithm and obtain results that are within a few per- 
cent of theoretical predictions. Figure 1 shows the tip 
velocity of the dendrite versus time for computations at 
undercoolings of A = 0.65 and 0.55. For all of these sim- 
ulations D — 1, d Q = 0.5, = 0, R = 15, and e = 0.05. 
For the A = 0.65 simulation, L = 200, Ax = 0.2, and the 
timestep is chosen to be At = 0.002. For the A = 0.55 
simulation, L = 800, Ax = 0.4, and At = 0.008. To en- 
sure grid convergence, Ax and At were rehned until the 
steady-state tip velocity did not vary by more than 2%. 

We also compare our level set results with simula- 
tions of dendritic growth performed using a phase-held 
model. Although calculations using phase-held models 
have been compared with a steady-state theory, there 
have been no comparisons made between time-dependent 
phase-field calculations and time-dependent solutions of 
the sharp-interface equations for multi-dimensional den- 
dritic growth. The phase-held model calculations pre- 
sented here were performed using a special adaptive mesh 
algorithm, as described in Ref. pl| . The tip velocity 
data from the phase-held model and level set method 
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FIG. 1. Time evolution of the tip velocity for simulations 
at A = 0.65 and 0.55. 

at A = 0.55 are in excellent agreement with each other 
(within 3%), as shown in Figure 1. Similar agreement 
is found in the dendritic shapes for these simulations, 
presented at time=9400 in Figure 2. These comparative 
results, combined with the recent demonstration of the 
equivalence of various phase-held models |2t|| , provide an 
excellent foundation for the validity of the phase-field ap- 
proach in simulating solidihcation microstructures. 
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FIG. 2. Comparison of dendritic shapes computed from the 
level set method and phase-field model for A = 0.55, shown 
at time=9400. 

The results presented here so far have used D$ = Dl, 
where D$ and Dl are the diffusivities in the solid and 
liquid phases, respectively. With our level set algorithm, 
we can also investigate the more general case where the 
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diffusivities are unequal. We performed additional simu- 
lations at A = 0.65 with D s = 0.75, 0.5, 0.25 and while 
keeping Dl = 1. The only available benchmark for the 
case of non-symmetric diffusion is the linearized solvabil- 
ity theory of Barbieri and Langer p2[ , which predicts 



P 2 V 



1 + D S /D L 



(P 2 V) DC 



(6) 



where p and V are the steady-state tip radius and veloc- 
ity, respectively. The values of p 2 V obtained from the 
level set simulations are compared with Eq. 6 in Figure 
3. The fit is surprisingly good (an error of about 13% 
at Ds/Dl — 0), considering Eq. 6 was obtained from a 
linearized theory in the limit of small undercoolings. 
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O Level Set Method A=0.65 



FIG. 3. p 2 V for different values of Ds / 'Dl- The circles are 
data from level set simulations at A = 0.65. The solid line is 
the theoretical prediction of Barbieri and Langer fitted to the 
data point at Ds/Dl = 1. 

In conclusion, the level set method should be consid- 
ered as a viable alternative to the use of phase-field mod- 
els. We have used a level set algorithm that can pro- 
duce accurate calculations of dendritic growth which can 
be compared favorably with solvability theory as well as 
time-dependent phase-field model simulations. The level 
set method can also handle discontinuous material prop- 
erties easily, which is currently very difficult with the 
phase-field approach. However, we should note that our 
implementation is not at all efficient. The practical ap- 
plication of this method to more realistic systems will 
require some sort of adaptive technique. In the future, 
we would like to use more computationally efficient im- 
plementations of this algorithm and apply these methods 
to problems in directional solidification, where the ability 
to simulate unequal diffusivities is of great interest. 
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